#include "VolumeMesh.h"
#include "VolumeGridData.h"
#include "MarchingCubes/MarchingCubes.h"
#include "rapidjson/writer.h"
#include "rapidjson/stringbuffer.h"
#include "OCMain/StringConverter.h"
#include "OCMain/zip/gzstream.h"

namespace OCPP
{
	namespace Volume
	{

		//_____________________________________________________________________________
		// Compute data
		void compute_data(MarchingCubes& mc)
			//-----------------------------------------------------------------------------
		{
			float x, y, z;
			float sx, sy, sz;
			float tx, ty, tz;

			float r, R;
			r = 1.85f;
			R = 4;

			sx = (float)mc.size_x() / 16;
			sy = (float)mc.size_y() / 16;
			sz = (float)mc.size_z() / 16;
			tx = (float)mc.size_x() / (2 * sx);
			ty = (float)mc.size_y() / (2 * sy) + 1.5f;
			tz = (float)mc.size_z() / (2 * sz);

			for (int k = 0; k < mc.size_z(); k++)
			{
				z = ((float)k) / sz - tz;

				for (int j = 0; j < mc.size_y(); j++)
				{
					y = ((float)j) / sy - ty;

					for (int i = 0; i < mc.size_x(); i++)
					{
						x = ((float)i) / sx - tx;
						mc.set_data((float)(
							// cushin
							z * z * x * x - z * z * z * z - 2 * z * x * x + 2 * z * z * z + x * x - z * z - (x * x - z) * (x * x - z) - y * y * y * y - 2 * x * x * y * y - y * y * z * z + 2 * y * y * z + y * y,
							// sphere
							//            ( (x-2)*(x-2) + (y-2)*(y-2) + (z-2)*(z-2) - 1 ) * ( (x+2)*(x+2) + (y-2)*(y-2) + (z-2)*(z-2) - 1 ) * ( (x-2)*(x-2) + (y+2)*(y+2) + (z-2)*(z-2) - 1 )) ,
							//  plane
							//            x+y+z -3,
							// cassini
							//            (x*x + y*y + z*z + 0.45f*0.45f)*(x*x + y*y + z*z + 0.45f*0.45f) - 16*0.45f*0.45f*(x*x + z*z) - 0.5f*0.5f ,
							// blooby
							//           x*x*x*x - 5*x*x+ y*y*y*y - 5*y*y + z*z*z*z - 5*z*z + 11.8 ),
							//  chair
							//            x*x+y*y+z*z-0.95f*25)*(x*x+y*y+z*z-0.95f*25)-0.8f*((z-5)*(z-5)-2*x*x)*((z+5)*(z+5)-2*y*y ,
							// cyclide
							//            ( x*x + y*y + z*z + b*b - d*d ) * ( x*x + y*y + z*z + b*b - d*d ) - 4 * ( ( a*x - c*d ) * ( a*x - c*d ) + b*b * y*y ),
							// 2 torus
							((x * x + y * y + z * z + R * R - r * r) * (x * x + y * y + z * z + R * R - r * r) - 4 * R * R * (x * x + y * y)) *
							((x * x + (y + R) * (y + R) + z * z + R * R - r * r) * (x * x + (y + R) * (y + R) + z * z + R * R - r * r) - 4 * R * R * ((y + R) * (y + R) + z * z))),
							// mc case
							// - 26.5298*(1-x)*(1-y)*(1-z) + 81.9199*x*(1-y)*(1-z) - 100.68*x*y*(1-z) + 3.5498*(1-x)*y*(1-z)
							// + 24.1201*(1-x)*(1-y)*  z   - 74.4702*x*(1-y)*  z   + 91.5298*x*y*  z  - 3.22998*(1-x)*y*  z  ),
							// Drip
							//          x*x + y*y - 0.5*( 0.995*z*z + 0.005 - z*z*z ) +0.0025 ),  // -0.0754+0.01, -0.0025 + 0.01, grid 40^3, [-1.5,1.5]

							i, j, k);
					}
				}
			}
		}
		//_____________________________________________________________________________

		VolumeMesh::VolumeMesh()
		{
			mcPtr = std::make_shared<MarchingCubes>();
		}

		bool VolumeMesh::convert(String path, FloatVector isoValues, int clipZ, const VolumeGridData& volumeGridData)
		{
			const auto& extent = volumeGridData.getExtent();
			int dimZ = volumeGridData.getExtent().dimZ - clipZ;
			int dimY = volumeGridData.getExtent().dimY;
			int dimX = volumeGridData.getExtent().dimX;

			MarchingCubes mc;

			FloatVector& vecFloat = (FloatVector&)volumeGridData._vecFloat;

			mc.set_ext_data((float*)&vecFloat[0]);

			for (auto isoValue : isoValues)
			{
				//float isoValue = 15.0f + i * 5.0f;
				mc.set_resolution(dimX, dimY, dimZ);
				mc.init_all();

				mc.run(isoValue);

				String outputName = str_format("%s_%.1f", path.c_str(), isoValue);
				//mc.writePLY((outputName + ".ply").c_str());

				rapidjson::StringBuffer buffer;
				rapidjson::Writer<rapidjson::StringBuffer> writer(buffer);
				writer.StartObject();
				writer.Key("numLon");
				writer.Int(dimX);
				writer.Key("numLat");
				writer.Int(dimY);
				writer.Key("numLevel");
				writer.Int(dimZ);

				writer.SetMaxDecimalPlaces(6);
				writer.Key("startLon");
				writer.Double(extent.xMin);
				writer.Key("endLon");
				writer.Double(extent.xMax);
				writer.Key("startLat");
				writer.Double(extent.yMin);
				writer.Key("endLat");
				writer.Double(extent.yMax);

				writer.Key("vertex");
				writer.SetMaxDecimalPlaces(3);
				writer.StartArray();

				for (size_t k = 0; k < mc.nverts(); k++)
				{
					auto v = mc.vertices()[k];
					writer.Double(v.x);
					writer.Double(v.y);
					writer.Double(v.z);
				}

				writer.EndArray();

				writer.Key("indices");
				writer.StartArray();
				for (size_t k = 0; k < mc.ntrigs(); k++)
				{
					auto v = mc.triangles()[k];
					writer.Int(v.v1);
					writer.Int(v.v2);
					writer.Int(v.v3);
				}
				writer.EndArray();

				writer.EndObject();

				bool gzipCompress = true;
				String outputFileName = outputName + ".json";
				if (gzipCompress)
				{
					outputFileName += ".gz";
					ogzstream ofs(outputFileName.c_str());
					ofs << buffer.GetString();
					ofs.close();
				}
				else
				{
					std::ofstream ofs(outputFileName);
					ofs << buffer.GetString();
					ofs.close();
				}

				mc.clean_all();
			}

			return true;
		}

		FloatVector VolumeMesh::run(float isoValue,const VolumeGridData& volumeGridData, const VolumeOptions& options)
		{
			int dimZ = volumeGridData.getExtent().dimZ;
			int dimY = volumeGridData.getExtent().dimY;
			int dimX = volumeGridData.getExtent().dimX;

			int multiple = 5;
			int bufferSize = 2046;
			size_t faceSize = bufferSize * bufferSize * 3 * multiple;
			unsigned int size = faceSize;
			float widthScale = dimX - 1;
			float heightScale = dimY - 1;
			float depthScale = dimZ - 1;
			if (depthScale == 0.0)
				depthScale = 1.0;

			FloatVector fv;
			fv.resize(size);

			mcPtr->set_ext_data((float*)&volumeGridData._vecFloat[0]);
			mcPtr->set_resolution(dimX, dimY, dimZ);
			mcPtr->init_all();
			mcPtr->run(isoValue);

			size_t offset = 0;
			size_t linmitOffset = offset + faceSize;
			for (size_t k = 0; k < mcPtr->ntrigs(); k++)
			{
				auto& tri = mcPtr->triangles()[k];
				auto& v0 = mcPtr->vertices()[tri.v1];
				auto& v1 = mcPtr->vertices()[tri.v2];
				auto& v2 = mcPtr->vertices()[tri.v3];
				if (offset + 9 >= linmitOffset)
					break;
				fv[offset++] = v0.x / widthScale;
				fv[offset++] = v0.y / heightScale;
				fv[offset++] = v0.z / depthScale;
				fv[offset++] = v1.x / widthScale;
				fv[offset++] = v1.y / heightScale;
				fv[offset++] = v1.z / depthScale;
				fv[offset++] = v2.x / widthScale;
				fv[offset++] = v2.y / heightScale;
				fv[offset++] = v2.z / depthScale;
			}

			mcPtr->clean_all();
			return fv;
		}

	}
}




